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摘要 田 谐 项 摄 动 是 分 析 法 轨道 预报 中 的 重要 部 分 , 其 中 包含 大 量 倾角 函数 及 其 偏 导数 的 计算 . 由 于 具有 精 
更 高 、 速 度 更 快 的 优点 , 倾角 函数 一 般 通 过 递 推 方法 计算 . 以 文献 中 提出 的 改进 Gooding 方 法 为 基础 , 将 其 给 
的 程序 稍 加 改进 , 在 计算 2-50 阶 倾角 函数 时 缩短 了 约 24% 的 计算 时 间 . 考虑 到 分 析 法 预报 过 程 中 轨道 平 
-一 化 很 小 , 以 泰勒 展开 式 计算 倾角 函数 , 可 极 大 提高 计算 速度 , 较 大 程度 地 减 小 分 析 法 预报 耗 时 ， 且 引 力 场 阶 次 
, 减 小 幅度 越 大 , 取 50 阶 时 预报 耗 时 缩短 了 48%. 另 一 方面 ,以 2 阶 展开 式 计算 倾角 函数 时 , 与 改进 Gooding 法 
HE, 分 析 法 预报 星 历 偏 差 很 小 . 对 于 500 km 高 度 的 低 轨 卫 星 , 分 别 以 改进 Gooding 法 和 2 阶 泰勒 展开 式 计算 贷 
(a ) 函数 , 预报 3 天 , 当地 球 引力 场 阶 次 不 高 于 50 时 , 二 者 预报 星 历 偏差 RMS (Root Mean Square) 低 于 1 mm, H. 
- 随 着 轨道 高 度 的 增加 , 预报 星 历 偏差 RMS 逐 渐 减 小. 
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2 分 析 法 轨道 预报 由 于 具有 较 大 的 速度 优势 而 ”2 惩 倾角 函数 的 3 个 指标 , 其 取 值 范围 为 2 < 1 < L, 
闪 被 广泛 应 用 于 近 地 天 体 动力 学 研究 中 , 在 构造 轨道 «m et op < 名 地 于 引力 了 次 .之 所 
Bl uu: TzsakPI A tH SE Jn E IA SCIL, (D), Hin. jA 
C o 形 引 力 、 日 月 第 三 体 引力 、 大 气 阻力 、 太 阳光 压 RO dn 
NE C Opa E 在 此 基础 上 Alanf 改 写 了 相应 符号 而 保持 公式 结 
V "bes á id 构 不 变 , 得 到 迄今 为 止 最 泛 用 的 倾角 函数 表达 式 
K, 与 二 体 引 力 相 比 量 级 约 为 DO(10-3). 地 球 非 球形 公式 如 下 ; 
引力 摄 动 函数 一 般 以 地 固 坐 标 系 中 卫星 的 球 坐 标 
进行 计算 , 若 将 球 坐 标 转换 为 卫星 轨道 根 数 ,那么 Fim) = 
函数 表达 式 中 可 分 离 出 只 包含 轨道 倾角 的 倾角 函 nu Di pe 2p ): 
数 , 而 倾角 函数 是 田 谐 项 摄 动 解 的 重要 组 成 因子 . p)! l-m-v 
在 近 地 天 体 动力 学 领域 , Kaulam 最 先 提 出 倾 CD (1) 
角 函 数 的 概念 并 记 为 wy( 了 ), 表达 式 为 关于 sin 了 和 2 
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其 中 v 为 求 和 指标 ， 

尽管 此 公式 已 足够 简洁 , 但 在 田 谐 项 摄 动 解 中 ， 
倾角 函数 需要 多 重 求 和 , 在 引力 场 阶 次 较 高 时 计算 
效率 低 , 且 会 出 现 损失 精度 的 情况 . Goodingl* 引 提 
出 另 一 种 递 推 计算 方法 , 将 归 一 化 倾角 函数 分 解 为 
两 个 函数 的 乘积 , 公式 如 下 : 


F5) = Vim) Anl), (2) 
其 中 Vs 是 关于 sin7 和 cos7 的 多 项 式 , ii AF. 则 是 关 
于 sin 志 和 cos 了 的 多 项 式 , 并 可 以 通过 三 项 递 推 公式 
计算 , Gooding jt ANu. FE E EE 关系 为 : 


(3) 


其 中 Ni 为 归 一 化 因子 , hk = 1 一 2p. 值得 注意 的 
是 倾角 函数 第 3 个 指标 无 论 使 用 p 还 是 k 都 可 以 ， 本 
质 上 没有 差别 , 而 Gooding 等 回 指出 虽然 & 取 值 不 连 
续 , 但 倾角 函数 关于 大 的 正 负 具有 对 称 性 , 在 递 推 计 
算 时 可 借 此 减少 计算 量 . 此 外 Gooding 等 回 对 递 推 
算法 做 出 进一步 优化 , 并 将 VS 分 离 出 不 含 倾角 的 
因子 WE 以 递 推 方法 计算 , 解决 了 在 高 阶 次 情况 下 
出 现 上 洪 和 下 溢 的 问题 . 同时 Gooding 等 回 还 给 出 
了 倾角 函数 1 阶 、2 阶 偏 导数 的 计算 公式 , 并 证 明了 
按照 (2) 式 计算 倾角 函数 时 效率 最 高 且 具 有 高 稳定 
因此 Gooding 递 推 公式 得 到 了 广泛 应 用 ， 成 为 
其 他 计算 方法 的 衡量 标准 . 

在 文献 [4-5] 中 Gooding 等 人 并 没有 详细 给 出 递 
推 公式 的 证 明 ,， 吴 连 大 等 6 分 别 利 用 超 几何 级 数 
的 递 推 关系 和 Jacobi 多 项 式 的 递 推 关系 , 完成 了 公 
式 推导 ,并 指出 其 实质 为 Jacobi 多 项 式 的 递 推 . 之 
后 吴 连 大 等 8- 在 Gooding 弟 推 方法 的 基础 上 提出 
计算 倾角 函数 的 改进 Gooding 方 法 , 按照 文中 给 出 
的 标准 递 推 过 程 计 算 44 ,而 对 于 Vi 则 将 Gooding 欠 
提出 的 表达 式 进 行 改 写 并 直接 计算 , 同时 避免 了 出 
现 上 洲 和 下 溢 的 情况 . 测试 结果 表明 改进 Good- 
ing 方 法 计算 速度 更 快 、 计 算 精 度 和 稳定 性 也 有 
所 提高 ， 男 一 方面 倾角 函数 的 计算 程序 也 更 加 简 
洁 ， 因 此 在 分 析 法 轨道 预报 过 程 中 可 以 使 用 改进 
Gooding 方 法 计算 倾角 函数 . 
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在 分 析 法 轨道 预报 过 程 中 地 球 非 球形 引力 摄 
动 解 的 计算 占据 了 较 大 比重 , 而 田 谐 项 摄 动 分 析 
解 中 需要 大 批量 计算 倾角 函数 ， 因 此 倾角 函数 的 
计算 速度 很 大 程度 上 影响 了 分 析 法 预报 耗 时 .本 
文 第 2 节 以 改进 Gooding 方 法 为 基础 , 对 附加 的 For- 
tran 程 序 进 行 改进 , 进一步 提升 了 大 批量 计算 倾角 
函数 的 效率 , 为 了 方便 , 本 文 称 此 方法 为 网 格 Good- 
ing 法 . 

一 般 情况 下 ， 对 地 球 卫 星 作 轨道 预报 时 预报 
弧 段 不 会 太 长 , 预报 过 程 中 轨道 倾角 的 变化 量 人 J 
很 小 ， 因 此 可 以 将 倾角 函数 在 初始 倾角 了 Hh 处 作 泰 
勒 展开 ， 以 1I 阶 或 2 阶 展 开 式 计算 倾角 函数 .倾角 
函数 展开 式 为 关于 A7 的 线性 或 2 次 多 项 式 ,与 改进 
Gooding 法 相 比 , 一 方面 计算 效率 得 到 大 幅 提升 , 但 
另 一 方面 展开 式 存在 一 定 的 截断 误差 , 从 而 导致 分 
析 法 预报 精度 受到 影响 . 如 果 展 开 式 误差 带 来 的 轨 
道 预报 结果 偏差 在 一 定 范围 内 , 那么 在 分 析 法 预报 
过 程 中 可 以 使 用 展开 式 计算 倾角 函数 , 从 而 很 大 程 
度 地 提高 预报 速度 ， 对 大 量 目标 作 轨 道 预报 时 利 
用 此 方法 可 极 大 减轻 工作 量 . 本 文 第 3 节 分 析 了 不 
同 初始 倾角 、 地 球 引 力 场 阶 次 、 倾 角 变 化 量 对 倾 
角 函 数 展开 式 误差 的 影响 , 比较 了 分 别 精确 至 1 阶 、 
2 阶 时 展开 式 的 计算 速度 以 及 精度 . 第 4 节 设 计 阁 干 
算 例 , 在 分 析 法 轨道 预报 中 分 别 使 用 网 格 Gooding 
法 、1 阶 、2 阶 展开 式 计算 倾角 函数 ， 以 改进 Good- 
ing 方 法 为 基准 , 分 析 预 报 星 历 偏差 以 及 预报 耗 时 ， 
其 中 预报 星 历 为 J2000 平 赤道 坐标 系 下 的 直角 坐标 . 


倾角 函数 递 推算 法 介绍 与 程序 优 
化 

由 于 吴 连 大 等 B79 提出 的 改进 Gooding 方 法 具 
有 计算 速度 快 、 计 算 精 度 和 稳定 性 高 、 程 序 简洁 
的 优点 , 故 可 以 直接 应 用 到 分 析 法 轨道 预报 中 , 而 
在 实际 应 用 中 发 现 弟 推 程序 有 进一步 提升 速度 的 
空间 . 根据 田 谐 项 摄 动 解 公式 可 知 预 报 过 程 中 需要 
计算 大 批量 倾角 函数 , 参考 刘 林 19 书 中 公式 : 


c(t) — y ` `> »» AOImpg ; 


122 m1 p—0 q——oo 


2 


(4) 
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其 中 c 包 的 表 示 轨 道 根 数 的 2 阶 短 周 期 项 ，Acuma 
中 包含 倾角 函数 mp 及 其 1 阶 偏 导数 本 ;,,， 刘 林 [9 
给 出 了 田 谐 项 摄 动 解 的 具体 公式 , 由 于 公式 较 长 , 
本 文 不 再 罗列 . 需要 注意 的 是 改进 Gooding 法 计算 
结果 为 归 一 化 倾角 函数 , 而 常用 的 分 析 解 中 倾角 函 
数 为 非 归 一 化 的 , 所 以 实际 计算 时 需 对 其 进行 转换 . 
根据 (4) 式 ,由 ! < 工 可 知 倾角 函数 计算 次 数 取 
决 于 地 球 引力 场 阶 次 , 递 推 次 数 N 可 由 下 式 计算 ， 


(5) 


 LQ-1042) , 
: 


这 表示 至 少 需 要 循环 调用 NN 次 递 推 程序 , 而 根据 改 
进 Gooding 法 原理 可 知 在 递 推 过 程 中 可 计算 出 若干 
倾角 函数 , 如 果 能 够 记录 下 这 些 值 , 那么 在 循环 过 
程 完全 可 以 跳 过 这 次 计算 , 从 而 减少 程序 调用 次 数 . 
当 引 力 场 阶 次 较 高 时 , 将 节省 不 少 倾角 函数 的 计算 
时 间 . 

同 Gooding 递 推 方法 一 样 ， 吴 连 大 等 &9 将 倾 
角 函 数 分 解 为 (2) 式 的 形式 ， 按 照 标准 过 程 递 推 计 
算 44 ,不 同 的 是 改写 了 V 表达 式 并 直接 计算 , 因 


此 在 批量 计算 时 只 需 考 虑 44# 的 程序 改进 , 其 递 推 
公式 为 : 
(J — 1)(J — m)(J + k) AS, (2) 
= (2j - 1)[J(J — 1)e — mk|AS (2) 
J(J4m-1)(J-k-1)45 44 (z), (6) 


其 中 x = cosI, 0 < J <l, 公式 沿 着 ! 递 推 , 递 推 初 
EN AS, FILAS ss 其 中 g = max(m, k). 该 递 推 关 
系 成 立 的 条 件 为 k > 0, 在 上 < 0 时 可 以 利用 倾角 函 
数 的 对 称 性 直接 计算 , 如 (7) 式 所 示 : 


Fa) ec) E D. 


(7) 


根据 (6) 式 可 知 ， 倾 角 函 数 在 沿 着 /方向 递 推 过 
程 中 可 得 到 中 间 值 4 (q < z < D), 由 (2) 式 和 (4) 式 


函数 展开 及 其 用 
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1 谐 项 摄 动 分 析 解 计算 中 所 


可 知 这 些 中 间 值 正 是 | 


需要 的 , 而 对 递 推 程序 的 改进 原理 为 将 其 记录 下 来 
以 避免 在 循环 过 程 中 重复 计算 . 


递 推 程序 改进 的 具体 方法 为 将 倾角 函数 以 3 维 
数组 形式 作为 函数 接口 . 在 吴 连 大 等 &-9] 给 出 的 程 
序 中 ,函数 输入 接口 为 轨道 倾角 和 3 个 指标 !/、m、 
k, 输出 接口 为 倾角 函数 及 其 偏 导 数 . 修改 之 后 输入 
接口 为 倾角 、 地 球 引 力 场 阶 次 L, 输出 接口 为 3 维 数 
组 形式 的 倾角 函数 及 其 偏 导 数 , 在 函数 内 部 A%, 同 
样 为 3 维 数 组 . 该 算法 的 原理 很 简单 ,规定 2 < 1 < L, 
1mcazl0zpcxl, 在 沿 着 :/、m、p 进 行 3 重 循 
环 过 程 中 将 p 转 换 为 k, 以 改进 Gooding 法 计算 倾角 
函数 , AR o FE 分 别 赋 给 3 维 数组 A[1][mj[p] 和 
[m][p]. SAGT SEAT, ZARTAN [m] [p EET ES 3, 
如 果 已 经 被 赋值 ， 则 跳 过 该 过 程 , 并 且 在 计算 倾角 
函数 时 直接 引用 AlN[mj[pl. 改进 之 后 仍 需 计算 和 NN 次 
VE, 而 4%, 的 计算 次 数 减 少 , 因此 倾角 函数 的 计算 
耗 时 有 所 缩减 , 且 随 着 阶 次 的 增加 , 耗 时 减少 比例 
逐渐 增 大 . 由 于 该 算法 本 质 为 使 用 改进 Gooding 法 
计算 倾角 函数 并 在 3 维 网 格 中 赋值 ， 因 此 精度 并 未 
改变 , 为 了 方便 , 本 文中 称 其 为 网 格 Gooding 法 . 

吴 连 大 等 gE9 测 试 了 1-200 阶 倾角 函数 的 计算 
时 间 , 与 Gooding 方 法 相 比 , 使 用 此 方法 时 效率 提 
升 了 41%. 本 文 设计 阁 干 算 例 , 对 改进 后 程序 的 计 
算 效 率 进 行 评估 ， 考 虑 到 在 分 析 法 预报 中 地 球 引 
力 场 取 50 阶 即 可 获得 高 精度 星 历 因此 选取 10、 
50, 每 组 设置 10 个 倾角 , EXX0.1*. 109. 209. 
90°, 分 别 按照 改进 Gooding 法 和 网 格 Gooding 法 批 
量 计算 倾角 函数 及 其 偏 导数 , 且 重 复 计 算 10 次 , 并 
记录 平均 耗 时 , 其 中 两 种 算法 都 以 C+ 十 语言 实现 . 
表 1 给 出 计算 机 配置 以 及 编译 设置 ， 所 有 算 例 均 
以 Release 模 式 进 行 测 试 . 相 比 之 下 , BELAK, 网 
格 Gooding 法 的 计算 效率 逐渐 提高 , L 取 10、50 时 计 
AFER] 2) Aki T 18.896. 23.896 ( 若 将 此 方法 应 用 
到 分 析 法 预报 中 , 预报 耗 时 可 参考 本 文 第 4 节 中 的 
表 5). 
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表 1 计算 机 配置 以 及 编译 设置 


Table 1 Computer configuration and compilation 


settings 


Test Environment Configuration 


CPU Intel (R)Core (TM)i7-8565U 
CPU Clock Speed 1.80 GHz 


Integrated Development . . 
Visual Studio 2019 


Environment 
Compiler LLVM-clang-cl 
Solution Configuration Release 
Optimization Level O2 


3 ”倾角 函数 的 展开 式 分 析 

3.1. 倾角 函数 展开 式 

根据 (1) 式 可 知 在 3 个 指标 不 变 时 倾角 函数 为 关 
于 sin Lflcos 4 的 多 项 式 , 若 已 知 初始 倾角 石 处 的 个 
角 函 数 及 其 偏 导数 , 则 可 以 在 五 附近 进行 泰勒 展开 ， 
以 展开 式 计算 倾角 函数 及 其 1 阶 偏 导数 ， 如 下 式 所 


ZN: 


FE (I) = FE (2o) + F,EQo)(I — Io)+ 
1 pk e 2 = 3 
g im (o) Io) + O(I1 — Io)”, (8) 
REO = FE o) + FzEo)U — Io)- 
OU — 19). 


通常 情况 下 ， 分 析 法 预报 的 弧 段 较 短 ， 轨 道 
平 倾角 变化 幅度 一 般 较 小 , 在 3 天 弧 段 以 内 基本 不 
超过 0.01*， 因 此 在 分 析 法 预报 中 能 够 对 倾角 函数 
作 泰 勒 展开 ,以 如 时 刻 的 倾角 函数 (7o) 及 其 1 阶 、 
2 阶 偏 导 数 为 初 值 , 在 t 时 刻 使 用 展开 式 计 算 倾 角 函 
数 . 倾角 函数 及 其 偏 导 数 初 值 通过 改进 Gooding 法 
计算 , 由 于 只 提供 2 阶 以 内 的 初 值 , 因此 倾角 函数 展 
开 式 最 高 可 以 精确 到 2 阶 , 而 偏 导数 只 能 展开 至 1 阶 . 

以 泰勒 展开 式 计算 倾角 通 数 可 大 幅 提 高 计算 
速度 , 从 而 有 效 缩短 分 析 法 预报 耗 时 . 与 改进 Good- 
ing 法 相 比 , 尽管 使 用 网 格 Gooding 法 计算 大 批量 倾 
角 函 数 能 在 保持 精度 不 变 的 同时 提高 计算 速度 , 但 


在 轨道 预报 过 程 中 倾角 函数 的 计算 只 是 其 中 一 部 
分 ,因此 预报 速度 虽 有 所 提升 , 但 并 不 明显 . 展开 
式 进 一 步 大 幅 减 小 了 倾角 函数 的 计算 量 , 经 测试 在 
地 球 引 力 场 取 50 阶 时 , 与 网 格 Gooding 法 对 比 , 1 阶 、 
2 阶 展开 式 的 计算 耗 时 均 缩短 了 99% 以 上 ， 因 此 不 
难 理解 以 展开 式 计算 倾角 函数 能 够 显著 提高 分 析 


法 预报 效率 . 
虽然 泰勒 展开 方法 可 以 有 效 提高 计算 效率 , 但 


是 计算 精度 却 受 到 一 定 影响 . 由 (8) 式 可 知 , 展开 式 
存在 截断 误差 , 倾角 函数 精度 降低 , 从 而 导致 预报 
结果 产生 偏差 . 只 有 当 截 断 误 差 较 小 、 预 报 星 历 偏 
差 不 足以 影响 预报 精度 时 , 才 可 以 将 倾角 函数 展开 
式 应 用 到 分 析 法 轨道 预报 中 . 


3.2 ”展开 式 误 差分 析 

吴 连 大 等 B79 已 经 验证 改进 Gooding 法 的 精度 
相 比 Gooding 法 更 高 一 些 , 而 网 格 Gooding 法 计算 
结果 与 改进 Gooding 法 一 致 , 因此 可 以 将 其 作为 计 
算 倾角 函数 的 精度 对 比 标准 . 倾角 函数 展开 式 误差 
按照 下 式 计算 : 


e1 = |F Qo + A1) — | Ff 49) + Fr I9)A1] 


, 


ea = | Fs Qo + A1) — | Fi (o) + FAC) A+ 


+ 


d= . 


(9) 


其 中 ej 为 倾角 函数 1 阶 展开 误差 , sz 为 2 阶 展 开 误差 ， 
2/ 为 倾角 函数 偏 导 数 1 阶 展开 误差 , 其 中 各 项 初 值 通 
过 网 格 Gooding 法 计算 得 到 . 
3.2.1. Hl 

根据 (9) 式 可 知 倾角 函数 展开 式 的 误差 主要 与 
初始 倾角 石 、 倾 角 变 化 量 AT 以 及 3 个 倾角 函数 指标 
有 关 , 为 了 较为 全 面 地 分 析 倾 角 函 数 展开 式 的 精度 
与 这 3 种 影响 因素 的 关系 , 我 们 设计 了 若干 算 例 如 
表 2 所 示 . 其 中 地 球 引 力 场 阶 次 取 为 50， 指 标 / 取 值 
范围 是 2 < 1 < 50, 初始 倾角 五 从 1" 起 每 隔 1" 取 值 


Fi (Io AT) — [Fs Co) + Ft Qo)A1] 


pu 
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至 179。， 此 外 也 设置 了 0.01? 及 179.99。， 总 共 181 组 
取 值 , 倾角 变化 量 和 7 分 别 取 为 0.1*>、0.01*、0.001°， 
倾角 函数 展开 式 精 确 至 1 阶 、2 阶 , 分 别 使 用 泰勒 展 
开 方 法 和 改进 Gooding 法 计算 倾角 函数 , 通过 这 些 
算 例 分 析 其 计算 结果 偏差 . 在 实际 应 用 场景 中 , 对 
于 近 圆 轨道 , 当地 球 引 力 场 阶 次 取 50, 使 用 完整 力 
模型 时 ， 以 分 析 方 法 预报 3 天 后 的 平 倾角 变化 量 基 
本 小 于 0.01*， 因 此 AT 选 择 0.01° 为 标准 . 倾角 函数 
有 3 个 指标 , 这 里 不 再 详细 讨论 m、p, 仅 计 算 各 
阶 次 ! 的 倾角 函数 误差 最 大 值 . 


表 2 倾角 函数 展开 式 误差 分 析 算 例 参数 设置 


Table 2 Parameters setting for the error analysis 


test of the inclination function expansion 


Parameters Values 
L 50 
l 2—50 
Io/ 0.01, 1-179, 179.99 
AI/ 0.1, 0.01, 0.001 
Expansion Order 1,2 


计算 结果 

由 于 倾角 函数 1 阶 、2 阶 展开 计算 速度 差别 很 
小 , 而 2 阶 展 开 式 精度 更 高 , 因此 主要 分 析 不 同 初始 
倾角 和 倾角 变化 量 对 2 阶 展 开 式 误差 的 影响 , Jf fig 
单 对 1 阶 、2 阶 展开 式 误差 进行 比较 . 2 阶 展 开 式 的 
误差 计算 结果 见 图 1, 分 别 表 示 A7 取 不 同 角度 时 侨 


3.2.2 


函数 展开 及 其 
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差分 布 与 上 述 规律 基本 一 致 , 只 是 误差 更 大 . 图 2 为 
倾角 函数 1 阶 展 开 误 差 结果 , 表示 A7 取 为 0.01" 时 展 
开 式 误差 最 大 值 的 等 高 线 图 , 误差 取 对 数 , 横 轴 为 
指标 l, 纵 轴 为 初始 倾角 J0. 根据 图 1 (d) 与 图 2 结果 ， 
倾角 函数 偏 导 数 的 1 阶 展开 误差 约 为 倾角 函数 1 阶 
展开 误差 的 30 信 . 

通过 图 1 中 6 个 子 图 的 对 比 可 知 , 倾角 函数 及 其 
导数 误差 随 着 AT 增 加 而 增 大 , 并 且 A7 每 增加 10 
Z, 倾角 函数 误差 增 大 约 1000 倍 , 而 偏 导 数 误差 增 
大 约 100 倍 ， 这 与 展开 式 的 截断 误差 是 AI 的 n 次 窜 
相关 , 为 便于 理解 , 将 (9) 式 整理 为 如 下 形式 : 


E: 


1 =t 
e1 = SFE (oA +OP), 


1 — HI 
e2 = S Bt (1)AT + O(AT?), 


1 — Ht 
gm 5 Fm (AT T O(AD). 


(10) 


据 此 可 知 倾角 函数 及 其 偏 导数 的 展开 式 误 差分 布 
和 2 阶 或 3 阶 偏 导数 分 布 是 一 致 的 . 

根据 计算 结果 显示 ， 倾 角 函 数 展开 误差 受 指 
标 ! 和 倾角 变化 量 AI 的 影响 较 大 , ! 越 大 误差 越 大 ， 
而 从 (10) 式 可 以 直观 看 出 误差 与 具有 窜 次 关系 . 
另外 初始 倾角 对 误差 影响 较 小 , 尤其 是 1 < 20 时 , 随 
着 了 0 增加 , 误差 几乎 没有 变化 . 通过 对 比 图 1 (c) 和 
图 2 可 知 倾角 函数 2 阶 展开 相 比 于 1 阶 展开 误差 减 小 
了 约 3 个 量 级 ， 上 述 分 析 全 部 基于 误差 最 大 值 进 
IT, 尽管 该 结果 不 能 完全 体现 出 倾角 函数 展开 误 


角 函 数 2 阶 展开 式 ( 左 列 ) 及 其 偏 导数 1 阶 展开 式 误差 
最 大 值 ( 右 列 ) 的 等 高 线 图 ， 其 中 误差 值 取 对 数 , 横 
为 指标 l， 纵 轴 为 初始 倾角 10. 由 图 1 可 知 随 着 ! 增 
D, 倾角 函数 2 阶 展开 误差 逐渐 增 大 ， 而 根据 等 高 


差 在 (l,m,) 空 间 的 分 布 特征 , 但 由 于 在 田 谐 项 摄 
动 分 析 解 计算 时 倾角 函数 误差 最 大 值 更 能 反映 对 
计算 精度 的 影响 , 所 以 上 述 结 果 已 经 足够 用 来 分 
析 倾 角 函 数 展开 式 的 适用 场景 了 . 对 于 预报 弧 段 较 


分 布 密度 可 知 误差 增 大 速度 逐渐 降低 . 另 一 方面 
EL < 20 时 ,误差 大 小 几乎 不 受 初始 倾角 的 影响 ， 
| > 20 时 ,误差 随 着 初始 倾角 增加 而 减 小 ， 当 超 
过 90° 之 后 又 逐渐 增 大 , 并 且 很 明显 可 以 看 出 等 高 
RRT Io = 90° 对 称 . 倾角 函数 偏 导 数 1 阶 展开 的 误 


È 
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短 、 平 倾角 变化 量 较 小 的 轨道 , 当地 球 引力 场 取 较 
低 的 阶 次 时 , 在 分 析 法 预报 中 可 以 尝试 将 倾角 函数 
以 2 阶 展开 式 计算 , 此 时 预报 误差 相对 较 小 , 不 过 实 
际 应 用 结果 需要 使 用 大 量 算 例 进行 测试 才能 进 
步 分 析 . 
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Maximum lge' 
b: AJ 7 0.1? 


lg(error) 


f: AI = 0.001? 


图 1 倾角 函数 2 阶 展开 式 及 其 


偏 导数 1 阶 


展开 式 的 误差 最 大 值 等 高 线 图 


Fig.1 Contour diagram for maximum error of 2nd-order expansion of the inclination function and 1st-order expansion of the 


partial derivative 
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Maximum lge 
AT = 0.01° 


lg(error) 


图 2 


函数 1 阶 展开 式 的 误差 最 大 值 等 高 线 图 


Fig.2 Contour diagram for maximum error of lst-order 


expansion of the inclination function 


4 ”分 析 法 轨道 预报 误差 分 析 

分 析 法 轨道 预报 中 使 用 展开 式 计算 倾角 函数 
可 有 效 提高 预报 速度 ,但 不 可 忽视 展开 式 存在 较 
大 误差 , 而 由 (4) 式 可 知 在 计算 田 谐 项 摄 动 分 析 解 
时 误差 会 在 三 重 求 和 过 程 中 逐渐 累积 ， 并 最 终 导 
致 预报 结果 产生 偏差 . 虽然 根据 图 1 与 图 2 可 以 初步 
分 析 倾 角 函 数 展开 式 的 适用 场景 , 但 为 了 进一步 分 
析 实 际 应 用 中 的 预报 精度 , 需要 使 用 大 量 算 例 作 轨 
道 预报 测试 . 一 般 情况 下 分 析 法 预报 误差 比较 大 ， 
预报 24 h 的 误差 为 数 十 米 至 数 百 米 不 等 , 而 倾角 函 
数 展开 导致 的 预报 误差 并 没有 这 么 大 , 如 果 直 接 将 
展开 后 的 预报 结果 与 数值 法 预报 结果 比较 , 则 可 能 
被 分 析 法 预报 本 身 的 误差 掩盖 . 由 于 展开 式 误差 不 
确定 正 负 , 实际 计算 中 可 能 出 现 展开 式 的 预报 精度 
比 完整 算法 更 高 的 “假象 ”, 因此 以 数值 法 预报 星 历 
作为 标准 无 法 准确 评估 展开 式 的 预报 精度 . 因为 整 
个 预报 过 程 中 只 有 倾角 函数 的 计算 方法 不 同 , 所 以 
直接 比较 展开 前 后 预报 星 历 的 偏差 即 可 . 本 章 以 改 
进 Gooding 法 为 基准 , 分 别 将 倾角 函数 展开 式 精确 
到 1 阶 和 2 阶 , 对 预报 星 历 偏差 进行 分 析 , 男 一 方面 ， 
为 定量 比较 倾角 函数 展开 前 后 分 析 法 预报 的 速度 ， 
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设计 了 若干 算 例 , 并 记录 预报 耗 时 . 


4.4 ”精度 分 析 算 例 

倾角 函数 展开 式 误 差 的 影响 因素 为 初始 倾 
角 I0、 倾 角 变 化 量 和 AT 和 阶 次 1, 而 阶 次 ! 的 取 值 范围 
与 地 球 引 力 场 阶 次 有 关 ， 因 此 选择 不 同 的 初始 轨 
道 , 分 析 初 始 倾角 bh 和 地 球 引 力 场 阶 次 L 对 预报 星 
历 偏差 的 影响 . 由 于 AT 是 在 轨道 预报 过 程 中 计算 
的 , 无 法 通过 给 定 初 值 分 析 , 所 以 不 再 考虑 此 因素 . 
谐 项 摄 动 分 析 解 公式 中 存在 因子 (上 ) 所 以 
容易 理解 轨道 半 长 径 o 越 大 倾角 函数 展开 误差 对 预 
I 影响 就 越 小 , 因此 主要 选取 低 轨 算 例 进行 
， 同 时 加 入 少量 中 高 轨 的 算 例 , 算 例 的 参数 设 
置 如 表 3 所 示 . 分 析 法 预报 初始 历 元 按照 简约 儒 略 
日 设 为 59000, 预报 长 度 为 3 d, 步 长 为 30 min, L4 
别 设 5 组 取 值 ， 地球 参考 半径 为 6378.137 km, 轨道 
高 度 h 通 过 将 轨道 半 长 径 a 与 地 球 参考 半径 相 减 得 
到 , 在 300 km 与 40000 km 之 间 取 值 , 共有 46 条 低 轨 
与 14 条 中 高 轨 , 轨道 偏心 率 均 设 为 0.001, 初始 倾角 
则 在 90? 以 内 每 隔 5 选取 了 19 组 值 ， 轨道 升 交 点 经 
度 、 近 地 点 角 和 平 近 点 角 均 设 为 50°. 


表 3 分 析 法 预报 精度 分 析 算 例 参 数 设 置 


Table 3 Parameters setting for test of calculation 


accuracy of analytical orbit prediction 


Parameters Values 
Initial Epoch (MJD) 59000 
Propagate Length/d 3 
Propagate Stepsize/s 1800 
L 4, 10, 20, 30, 40, 50 
Earth Reference Radius/km 6378.137 


300-2000 (46); 


Height /km 
3000-40000 (14) 
Eccentricity 0.001 
Io/? 0.01, 5, 10, --., 90 
Ascending Node/? 50 
Argument of Perigee/? 50 
Mean Anomaly/? 50 


分 析 法 预报 所 用 摄 动 模型 为 完整 力 模型 , 包括 
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带 谐 项 、 田 谐 项 、 日 月 引力 摄 动 、 大 气 阻 力 摄 动 
和 光 压 摄 动 等 . 一 般 情况 下 为 了 发 挥 分 析 法 预报 
的 速度 优势 , 步 长 设置 为 与 预报 时 长 相等 , 而 本 次 
测试 将 步 长 设 为 30 min, 一 方面 是 为 了 得 到 较 密集 
的 星 历 , 便于 分 析 星 历 偏差 另 一 方面 可 以 更 频繁 
地 更 新 参考 星 历 ， 从 而 一 定 程度 上 提高 预报 精度 
此 外 轨道 预报 是 精密 定 轨 过 程 中 不 可 或 缺 的 一 部 
分 ， 而 精密 定 轨 中 观测 数据 的 历 元 间隔 也 决定 了 
步 长 是 较 小 的 . 由 于 倾角 函数 展开 误差 与 偏心 率 
无 关 , 因此 简单 地 选取 近 圆 轨道 进行 测试 , 偏心 率 
取 0.001. 根据 倾角 函数 展开 式 误 差分 析 结 果 可 知 
误差 关于 万 = 90? 对 称 , 因此 不 再 考虑 初始 倾角 大 
于 90° 的 算 例 . 

在 计算 田 谐 项 摄 动 解 时 可 能 存在 共振 现象 ， 
即 (4) 式 中 某 些 项 具有 小 分 母 问 题 (分 母 接近 于 0)， 
这 主要 取决 于 求 和 指标 !/、m、p、9g 和 卫星 平 运动 
速率 , 具体 表达 式 可 见 文献 [10]. 一 般 情况 下 地 球 引 
力 场 阶 次 越 高 , 满足 小 分 母 关 系 的 求 和 指标 越 多 ， 


民 项 也 就 越 多 , 另外 低 轨 卫 星 会 更 频繁 地 遇见 共 
振 现象 . 在 本 文 所 使 用 的 分 析 法 预报 算法 中 , 计算 


田 谐 项 时 所 有 满足 小 分 母 关 系 的 求 和 指标 不 再 参 
与 短 周期 项 的 计算 


测试 中 分 别 使 用 网 格 Gooding 法 、1 阶 展开 和 2 
阶 展开 这 3 种 方法 计算 倾角 函数 及 其 偏 导数 ,遍历 
以 上 轨道 预报 算 例 , 得 到 3 组 预报 星 历 , 即 各 时 刻 的 
位 置 矢量 : re(t), ri(t), ra(t), 接着 按照 下 式 分 别 计 
算 1 阶 展开 和 2 阶 展开 的 预报 星 历 偏差 RMS: 

?()- r0 


2 


其 中 mt 为 星 历 点 数量 ,根据 RMS 结 果 对 倾角 函数 
展开 式 在 轨道 预报 中 的 适用 性 进行 分 析 . 


4.2 ” 算 例 结果 分 析 

为 研究 地 球 引 力 场 阶 次 和 初始 倾角 与 预报 精 
度 的 关系 , 挑选 出 相同 轨道 高 度 的 算 例 , 对 预报 星 
历 偏差 RMS 作 图 , 并 为 了 直观 分 析 , 将 RMS 取 对 数 . 
另外 挑选 出 若干 初始 倾角 相同 的 算 例 , 比较 不 同 轨 
道 高 度 的 预报 星 历 偏差 RMS, 从 而 分 析 在 不 同 轨道 
高 度 下 , 倾角 函数 展开 对 预报 结果 的 影响 . 由 于 倾 
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角 函 数 2 阶 展开 精度 更 高 , 且 速 度 几乎 相同 , 因此 主 
要 对 2 阶 展开 的 预报 结果 进行 分 析 . 

对 于 轨道 高 度 在 2000 km 以 内 的 低 轨 算 例 , 一 
般 情况 下 地 球 引 力 场 阶 次 越 高 , 预报 星 历 偏差 RMS 
BRK, 只 是 增加 幅度 很 小 . L < 50 时 , 预报 星 历 偏差 
RMS 整 体 上 随 着 初始 倾角 的 增加 有 逐渐 减 小 的 趋 
势 , 这 和 倾角 函数 误差 的 变化 规律 相符 , 而 在 
90° 时 RMS 小 幅 增 大 . 以 500 km 高 度 的 轨道 为 例 , 倾 
角 函 数 2 阶 展开 的 预报 星 历 偏差 如 图 3 所 示 , 横 坐 标 
为 初始 倾角 石 ， 纵 坐标 为 预报 星 历 偏差 RMS, 并 取 
对 数 ， 图 例 表示 地 球 引 力 场 阶 次 L， 由 图 3 可 知 地 
球 引 力 场 在 50 阶 以 内 时 RMS 小 于 0.5 mm, 在 
85° 处 RMS 为 极 小 值 . 经 过 进一步 计算 和 分 析 发 现 ， 
由 于 各 阶 次 的 引力 场 系数 和 倾角 函数 表达 式 不 同 ， 
不 同 阶 次 的 RMS 极 小 点 存在 差异 , 例如 L = 3 时 极 
小 点 在 万 = 83.5° 左 右 , 因为 田 谐 项 分 析 解 为 各 阶 
次 计算 结果 的 和 , 所 以 图 3 中 极 小 点 为 综合 结果 , 并 
不 严格 在 90° 处 . 

对 于 中 高 轨 算 例 , 地 球 引力 场 阶 次 对 预报 误差 
影响 很 小 , 而 随 着 初始 倾角 增加 , 轨道 高 度 不 同时 ， 
预报 星 历 偏差 RMS 具 有 不 同 的 变化 规律 , 但 RMS 
整体 很 小 ， 比 如 轨道 高 度 为 10000 km 算 例 的 RMS 
随 初始 倾角 增加 而 逐渐 增 大 , 而 高 度 为 36000 km 的 
同步 轨道 算 例 的 RMS 却 有 所 浮动 , 如 图 4 所 示 , 不 过 
二 者 RMS 均 不 超过 10-7 m. 

另 一 方面 为 更 清晰 地 分 析 预 报 误 差 与 轨道 高 
度 h 的 关系 , WL = 50, 将 预报 星 历 偏差 RMS 关 
TER]. 轨道 高 度 小 于 和 大 于 2000 km, 结果 分 别 
如 图 5、 图 6 所 示 , 其 中 横 坐 标 为 h, 纵 坐 标 为 预报 星 
历 偏差 RMS 的 对 数值 ， 图 例 表示 初始 倾角 0. 可 以 
发 现 当 轨 道 高 度 小 于 2000 km 时 , RMS 随 轨道 高 度 
增加 而 逐渐 减 小 , 但 随 着 初始 倾角 增加 而 减 小 , 且 
均 小 于 1 mm. 而 当 轨 道 高 度 大 于 2000 km 时 , RMS 
先 随 轨道 高 度 增加 而 减 小 ， 然 后 逐渐 增 大 并 维持 
在 10-8 m 至 约 10-7 m 之 间 . 

通过 上 述 算 例 结果 可 知 , 在 分 析 法 轨道 预报 中 ， 
当地 球 引 力 场 不 高 于 50 阶 时 ， 以 2 阶 泰勒 展开 式 计 
算 倾角 函数 ， 与 改进 Gooding 法 相 比 , 低 轨 卫 星 的 
预报 星 历 偏 差 RMS 小 于 1 mm, 而 中 高 轨 卫 星 则 更 
小 ， 以 地 球 同 步 轨道 卫星 为 例 RMS 约 为 10-7 m. 
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在 实际 应 用 场景 中 , 分 析 法 预报 一 般 不 需要 使 用 高 
阶 次 的 地 球 引 力 场 模型 , 50 阶 已 经 足够 满足 精度 要 
K, 由 于 分 析 法 预报 24 h 的 位 置 误差 一 般 大 于 10 m 
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量 级 ,因此 可 以 认为 倾角 函数 2 阶 展开 不 会 影响 分 


析 法 预报 精度 . 
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图 3 500 km 高 度 轨道 在 取 不 同 初始 倾 


9 时, di 


函数 2 阶 上 万 


开 的 分 析 法 预报 星 历 偏差 


Fig.3 The ephemeris deviations of analytical orbit prediction for 500 km altitude orbit that have different initial inclination 


with inclination function calculated by 2nd-order expansions 
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36000km: 2nd-order expansion 
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图 4 36000 km 高 度 轨道 在 取 不 同 初 始 倾角 时 , (1 
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的 分 析 法 预报 星 历 偏差 


Fig.4 The ephemeris deviations of analytical orbit prediction for 36000 km altitude orbit that have different initial inclination 


with inclination function calculated by 2nd-order expansions 
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图 5 不 同 高 度 的 低 轨 倾角 函数 2 阶 展开 的 分 析 法 预报 星 历 偏差 


Fig.5 The ephemeris deviations of analytical orbit prediction for Low Earth Orbit at different altitude with inclination function 


calculated by 2nd-order expansions 
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L=50: 2nd-order expansion 
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图 6 不 同 高 度 的 中 高 轨 倾 角 函数 2 阶 展开 的 分 析 法 预报 星 历 偏差 


Fig.6 The ephemeris deviations of analytical orbit prediction for Medium-High Earth Orbit at different altitude with 


inclination function calculated by 2nd-order expansions 


4.3 ”分 析 法 预报 速度 分 析 取决 于 地 球 引 力 场 阶 次 的 大 小 以 及 预报 节点 的 数 

为 定量 比较 倾角 函数 展开 前 后 分 析 法 预报 的 。” 量 , 因此 轨道 半 长 径 不 需 考虑 , 引力 场 阶 次 分 别 取 6 
效率 , 设计 若干 算 例 进行 分 析 , 参数 设置 如 表 4 所 示 . ^r, 预报 时 长 和 预报 步 长 分 别 设 为 3 d、30 min, 共 
由 于 分 析 法 预报 耗 时 受 轨 道 半 长 径 影 响 很 小 , 主要 ”有 144 个 预报 节点 . 摄 动 模型 依然 为 完整 力 模型 , 轨 


Im 
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CI 
x 


道 偏心 率 为 0.001, 初始 倾角 在 1%-90° 之 间 间 隔 1° 取 耗 时 最 多 可 减少 约 50%. 


值 . "-—' 
测试 环 HAUKI, 以 Release 模 式 编译 程序 . Table 4 iu e. Mine cc NR 

分 别 使 用 改进 Gooding 法 、 网 格 Gooding 法 、1 阶 展 speed of analytical orbit prediction 

开 和 2 阶 展开 4 种 方法 计算 倾角 函数 , 在 工 取 不 同 值 Parameters Values 

时 , 分 别 对 不 同 初始 倾角 的 90 组 算 例 进行 预报 , id Initial Epoch (MJD) 59000 

录 预 报 耗 时 并 求 和 , 并 以 改进 Gooding 法 为 基准 , 计 Propagate Length/ d 3 

算 其 余 3 种 方法 的 耗 时 减 小 比例 . 统计 结果 如 表 5 所 Propagate Štepsizē/s 1800 

示 , 网 格 Gooding 法 使 得 分 析 法 轨道 预报 速度 有 小 L 4, 10, 20, 30, 40, 50 

EEF, 在 五 较 高 时 预报 耗 时 减少 13% 左 右 , 而 1 阶 、 batis Adele — 

2 阶 展开 对 预报 速度 的 提升 效果 十 分 显著 ， 且 二 者 Éciiidiy —— 

差别 很 小 , PERLA TRE EE IST J^] EGLI TE 03 or ii 

K, L = 50 时 达到 48% 左 右 . 因为 分 析 法 预报 中 除 | . 

了 倾角 函数 还 有 其 他 各 种 计算 , 因此 预报 耗 时 的 减 S MI i 

小 存在 上 限 , 根据 测试 结果 分 析 , 倾角 函数 的 计算 B pps dd 50 

耗 时 约 占 预报 总 耗 时 的 50%, 因此 当 工 增加 时 , 预报 ss 


表 5 使 用 不 同方 法 计算 倾角 函数 时 分 析 法 预报 耗 时 对 比 


Table 5 The comparison of run times of analytical orbit predictions with inclination function calculated 
by different methods 


Modified Gooding's X Gridding Gooding's 1st-order 2nd-order Time Consumption 
Method/s Method/s Expansion/s | Expansion/s Reduction Ratio 
4 1.501 1.593 1.176 1.174 —696 2296 2296 
10 7.5 1.183 5.584 5.592 4% 2696 2596 
20 47.041 42.042 30.319 29.675 11% 36% 37% 
30 151.762 132.276 88.268 89.483 13% 42% 41% 
40 359.814 315.556 197.973 196.378 12% 45% 45% 
50 719.97 622.329 370.176 371.64 14% 49% 48% 
5 结论 展开 式 应 用 到 分 析 法 轨道 预报 中 的 可 行 性 , 设计 了 


本 文 基于 改进 Gooding 法 调整 了 倾角 函数 的 计 大 量 算 例 , 在 预报 过 程 中 分 别 使 用 改进 Gooding 法 
SURUT, 使 其 能 够 一 次 性 作 大 批量 计算 , 在 计算 50 和 1 阶 、2 阶 展开 式 计算 倾角 函数 , 分 析 项 报 星 历 全 
阶 以 内 的 倾角 函数 时 ， 相 较 于 循环 调用 程序 , 使 用 ARMS. 根据 测试 结果 , 在 地 球 引 力 场 取 50 阶 以 内 
网 格 Gooding 法 耗 时 缩短 了 23.8%. 时 , 倾角 函数 2 阶 展开 后 ,， 低 轨 卫 星 的 预报 星 历 偏 

考虑 到 倾角 函数 计算 公式 为 多 项 式 , 可 以 进行 。 差 RMS 小 于 1 mm, 而 随 着 轨道 高 度 增加 ,RMS 旦 
泰勒 展开 , 而 在 短 弧 轨 道 预报 中 平 倾角 变化 量 一 般 ”” 减 小 的 趋势 . 一 般 在 分 析 法 预报 中 , 地 球 引力 场 阶 
很 小 , 所 以 本 文 对 倾角 函数 展开 式 的 计算 精度 进行 。 ”次 取 值 较 低 , 基本 不 会 超过 50 阶 , 因此 将 倾角 函数 
分 析 , 发 现 2 阶 展开 式 精度 更 高 . 为 了 研究 倾角 函数 作 2 阶 展开 并 不 会 影响 分 析 法 预报 精度 . 
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Expansion of Inclination Functions and Its Application in 
Analytical Orbit Prediction 
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(2 Key Laboratory of Modern Astronomy and Astrophysics, Ministry of Education, Nanjing 210023) 
(3 Institute of Space Environment and Aerospace Dynamics, Nanjing University, Nanjing 210023) 


ABsTRAcT Tesseral perturbation is a major part in analytical orbit propagation, which involves a sig- 
nificant amount of calculations of inclination functions and their partial derivatives. Recursion formula is 
widely applied for accuracy and efficiency purposes. Based on the modified Gooding's method proposed in 
the literature, an optimized program is proved to reduce 2496 of calculation time for all inclination func- 
tions up to degree/order 50 (L < 50). Furthermore, considering that the mean inclination varies slightly 
in orbit prediction for 1—3 days, the functions are expanded as Taylor series, which speeds up calculations 
significantly, and the time consumed during orbit prediction is reduced by 4896 while the tesseral pertur- 
bation is calculated up to degree/order 50. Moreover, with expansion of inclination functions up to the 
second order, the predicted ephemerides deviate slightly from those calculated using full recursion for- 
mula. For a 500 km low Earth satellite, propagating its orbits for 3 days using modified Gooding's method 
and second-order Taylor expansion respectively to calculate the inclination functions within degree/order 
50, the ephemeris deviation RMS (Root Mean Square) is less than 1 mm and decreases as the altitude 
increases. 


Key words celestial mechanics: inclination function, modified Gooding's recursion, Taylor expansion, 
analytical orbit prediction 
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